# This script produces Figures 4, 5, and 13 through 15 # 
# NOTE THAT YOU MUST INSTALL AN OLDER VERSION OF PANELMATCH
# USING PanelMatch.zip before running this script

rm(list=ls())
options(digits=4)
pkg <- c("PanelMatch", "ggplot2", "plm", 
         "matrixStats",
         "foreign", 
         "stringi", "stringdist",
         "DataCombine","lmtest", "multiwayvcov",
         "data.table",
         "dplyr", 
         "stargazer",
         "clubSandwich", "mediation",
         "doBy")
lapply(pkg, require, character.only = TRUE)

# load datasets
proprietary <- read.csv("proprietary.csv")
load("non_p.RData")

# merge 
d_sub1 <- merge(d_sub1, proprietary, by = c("county_id", "time"), all.x = T)

# creating the 5th and 6th LDVs
d_sub1 <- slide(data = d_sub1, Var = "log_transaction_area", GroupVar = "county_id",
                TimeVar = "time", slideBy = -5, NewVar = "log_transaction_area_l5")

d_sub1 <- slide(data = d_sub1, Var = "log_transaction_area", GroupVar = "county_id",
                TimeVar = "time", slideBy = -6, NewVar = "log_transaction_area_l6")

## Figure 4 ##
matches_M50 <- PanelMatch(formula = log_transaction_area ~ 
                            tour + 
                            log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                            log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                            log_expd_l1 + log_middle_school_l1 + 
                            msec2currentsec_l1 + 
                            msec2currentgvn_l1 + 
                            mayor2currentsec_l1 + 
                            mayor2currentgvn_l1 + 
                            msec2currentsec2_l1_missing, 
                          lag = 6, max.lead = 4,
                          unit.id = "county_id",
                          time.id = "time",
                          treatment = "tour",
                          method = "Maha",
                          M = 50,
                          weighting = FALSE,
                          qoi = "ate",
                          data = d_sub1)    

matches_M20 <- PanelMatch(formula = log_transaction_area ~ 
                            tour + 
                            log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                            log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                            log_expd_l1 + log_middle_school_l1 + 
                            msec2currentsec_l1 + 
                            msec2currentgvn_l1 + 
                            mayor2currentsec_l1 + 
                            mayor2currentgvn_l1 + 
                            msec2currentsec2_l1_missing, 
                          lag = 6, max.lead = 4,
                          unit.id = "county_id",
                          time.id = "time",
                          treatment = "tour",
                          method = "Maha",
                          M = 20,
                          weighting = FALSE,
                          qoi = "ate",
                          data = d_sub1)    

### ATT ###
covariates <- c("log_gdppc_l1", "log_gdp_l1", "log_fia_l1",
                "log_loan_l1", "log_hos_l1", "log_industrial_l1",
                "log_rev_l1", "log_expd_l1", "log_middle_school_l1",
                "msec2currentsec_l1", "msec2currentgvn_l1", 
                "mayor2currentsec_l1", "mayor2currentgvn_l1")

balance_DV_nonref_M50 <- PanelBCheck(matched_sets = matches_M50, covariate = NULL,
                                     qoi = "att", refinement = FALSE)$Mean_Balance$balance

balance_DV_ref_M50 <- PanelBCheck(matched_sets = matches_M50, covariate = NULL,
                                  qoi = "att", refinement = TRUE)$Mean_Balance$balance

balance_DV_ref_M20 <- PanelBCheck(matched_sets = matches_M20, covariate = NULL,
                                  qoi = "att", refinement = TRUE)$Mean_Balance$balance


balance_DVs <- list(balance_DV_nonref_M50, 
                    balance_DV_ref_M50,
                    balance_DV_ref_M20)

balance_DVs_pre <- list(balance_DV_nonref_M50[1:6], 
                        balance_DV_ref_M50[1:6],
                        balance_DV_ref_M20[1:6])

pdf(file = "output/Figure4.pdf", 
    width = 12, height = 7)
par(mfrow = c(1,3), 
    mar = c(2, 2, 2 , 1), oma = c(3, 5,4, 1))
ylim = c(-1, 1)
for (i in 1:length(balance_DVs_pre)) {
 
  plot(1:length(balance_DVs_pre[[i]]), 
       balance_DVs_pre[[i]], 
       ylim = ylim, 
       # ylab = ylab,
       # main = main,
       type = "l", xaxt = "n", yaxt = "n")     
  
  axis(side = 1, at = 1:length(balance_DVs_pre[[i]]), 
       labels = names(balance_DVs_pre[[i]]))
  axis(side = 2, at = seq(-1, 1, .25), 
       labels = seq(-1, 1, .25))
  abline(h = .25, lty = 5)
  abline(h = -.25, lty = 5)
  abline(h = 0, lty = 3)
  abline(v = 6, lty = 3)
}

mtext("Raw Data after Matching", side = 3, line = .6, at = 0.16636,
      outer = TRUE, cex = 1.2)

mtext("Refinement with \n Fifty Closest Control Units", side = 3, 
      line = .6, at = 0.502,
      outer = TRUE, cex = 1.2)

mtext("Refinement with \n Twenty Closest Control Units", side = 3, 
      line = .6, at = 0.835,
      outer = TRUE, cex = 1.2)

mtext("Standardized Mean Differences \n between Treated and Untreated", 
      side = 2, line = 0.4, at = 0.5,
      outer = TRUE, cex = 1.2)

mtext("Months Relative to the Administration of Treatment", side = 1, line = 1.385, at = 0.53,
      outer = TRUE, cex = 1.2)

dev.off()

##########################
# Figures 5 and 13 through 15
cores <- 4
## L = 4, Maha ##
Maha_M50_ATE <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 50,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M50_ATE, file = "output/Maha_M50_ATE.RData")

##
Maha_M20_ATE <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 

                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 20,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M20_ATE, file = "output/Maha_M20_ATE.RData")

##
Maha_M50_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 50,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M50_ATE_L6, file = "output/Maha_M50_ATE_L6.RData")

##
Maha_M20_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 20,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M20_ATE_L6, file = "output/Maha_M20_ATE_L6.RData")

#### Maha Figures M = 20 and 50 #####
load("output/Maha_M20_ATE.RData")
load("output/Maha_M50_ATE.RData")
load("output/Maha_M20_ATE_L6.RData")
load("output/Maha_M50_ATE_L6.RData")

all_results <- list(Maha_M50_ATE, Maha_M20_ATE, Maha_M50_ATE_L6, 
                    Maha_M20_ATE_L6)


pdf(file = "output/Figure5.pdf", 
    width = 10, height = 8)
par(mfrow = c(2,2), 
    mar = c(2, 3, .5 , 1), oma = c(3, 3,4, 1))

for (i in 1:length(all_results)) {
  total_periods <- length(all_results[[i]][[length(all_results[[i]])]]$lead)
  P <- rep(NA, total_periods)
  L_90 <- rep(NA, total_periods)
  U_90 <- rep(NA, total_periods)
  L_95 <- rep(NA, total_periods)
  U_95 <- rep(NA, total_periods)
  CI_90 <- .9
  CI_95 <- .95
  if (i == 1|i == 3){
    ylab = "Estimated Effect of Inspection"
  } else {
    ylab = ""
  }
  for (j in 1:length(all_results[[i]])) {
    results <- all_results[[i]][[j]]
    if (j == 1){
      P[1:length(results$lead)] <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      L_90[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_90
      U_90[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      L_95[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_95
      U_95[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
    } else {
      tmp <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      P[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      U_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_95[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
      U_95[results$lag + j] <- tmp[length(tmp)]
    }   
  }
  grey_gov <- length(P) - 5 - .5       
  xlim_gov <- length(P) - .5
  P[length(P)-5] <- NA
  L_90[length(L_90)-5] <- NA
  U_90[length(U_90)-5] <- NA
  L_95[length(L_95)-5] <- NA
  U_95[length(U_95)-5] <- NA
  ylim <- c(min(L_95, na.rm = T) * 1.1, 
            max(U_95, na.rm = T) * 1.1) 
  plot(0:(length(as.numeric(P))-1), P, xlab = "", ylab = ylab, 
       xaxt = "n",
       # yaxt = "n",
       xlim = c(-.5,xlim_gov),
       
       panel.first = rect(-10, -10, grey_gov, 10, border = NA,
                          col = "grey88"), 
       ylim = ylim, 
       pch = 16, cex = 1.4) 
  #  }
  
  # #  if (k == 1) {
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_95[k+1], U_95[k+1]))
  }
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_90[k+1], U_90[k+1]), lwd = 5,
          col = "grey")
  }
  points(0:(length(as.numeric(P))-1), P, 
         pch = 16, cex = 1.4)
  abline(h = 0, lty = 3)
  axis(side = 1, at = 0:(length(P)-1), 
       labels = -(length(P)-5):4)
}

mtext("Months Relative to the Administration of Treatment", side = 1, line = 1.385, at = 0.53,
      outer = TRUE, cex = 1.2)

mtext("Twenty Closest Control Units", side = 3, line = 1.385, at = 0.268,
      outer = TRUE, cex = 1.2)

mtext("Fifty Closest Control Units", side = 3, line = 1.385, at = 0.78,
      outer = TRUE, cex = 1.2)

mtext("Four Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.76,
      outer = TRUE, cex = 1.2)

mtext("Six Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.262,
      outer = TRUE, cex = 1.2)

dev.off()


## Pscore M = 50 , M = 20 ##
cores <- 4
Pscore_M50_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 50,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M50_ATE_L6, file = "output/Pscore_M50_ATE_L6.RData")

Pscore_M20_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 20,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M20_ATE_L6, file = "output/Pscore_M20_ATE_L6.RData")

Pscore_M50_ATE_L4 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 50,
                        weighting = FALSE,

                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M50_ATE_L4, file = "output/Pscore_M50_ATE_L4.RData")

Pscore_M20_ATE_L4 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 20,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M20_ATE_L4, file = "output/Pscore_M20_ATE_L4.RData")


#### Pscore Figures ####
load("output/Pscore_M20_ATE_L4.RData")
load("output/Pscore_M50_ATE_L6.RData")
load("output/Pscore_M20_ATE_L6.RData")
load("output/Pscore_M50_ATE_L6.RData")

all_results <- list(Pscore_M50_ATE_L6, Pscore_M20_ATE_L4, Pscore_M50_ATE_L6, 
                    Pscore_M20_ATE_L6)

pdf(file = "output/Figure14.pdf", 
    width = 10, height = 8)
par(mfrow = c(2,2), 
    mar = c(2, 3, .5 , 1), oma = c(3, 3,4, 1))

for (i in 1:length(all_results)) {
  total_periods <- length(all_results[[i]][[length(all_results[[i]])]]$lead)
  P <- rep(NA, total_periods)
  L_90 <- rep(NA, total_periods)
  U_90 <- rep(NA, total_periods)
  L_95 <- rep(NA, total_periods)
  U_95 <- rep(NA, total_periods)
  CI_90 <- .9
  CI_95 <- .95
  if (i == 1|i == 3){
    ylab = "Estimated Effect of Inspection"
  } else {
    ylab = ""
  }
  for (j in 1:length(all_results[[i]])) {
    results <- all_results[[i]][[j]]
    if (j == 1){
      P[1:length(results$lead)] <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      L_90[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_90
      U_90[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      L_95[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_95
      U_95[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
    } else {
      tmp <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      P[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      U_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_95[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
      U_95[results$lag + j] <- tmp[length(tmp)]
    }   
  }
  grey_gov <- length(P) - 5 - .5       
  xlim_gov <- length(P) - .5
  P[length(P)-5] <- NA
  L_90[length(L_90)-5] <- NA
  U_90[length(U_90)-5] <- NA
  L_95[length(L_95)-5] <- NA
  U_95[length(U_95)-5] <- NA
  ylim <- c(min(L_95, na.rm = T) * 1.1, 
            max(U_95, na.rm = T) * 1.1) 
  plot(0:(length(as.numeric(P))-1), P, xlab = "", ylab = ylab, 
       xaxt = "n",
       # yaxt = "n",
       xlim = c(-.5,xlim_gov),
       
       panel.first = rect(-10, -10, grey_gov, 10, border = NA,
                          col = "grey88"), 
       ylim = ylim, 
       pch = 16, cex = 1.4) 
  #  }
  
  # #  if (k == 1) {
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_95[k+1], U_95[k+1]))
  }
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_90[k+1], U_90[k+1]), lwd = 5,
          col = "grey")
  }
  points(0:(length(as.numeric(P))-1), P, 
         pch = 16, cex = 1.4)
  abline(h = 0, lty = 3)
  axis(side = 1, at = 0:(length(P)-1), 
       labels = -(length(P)-5):4)
}

mtext("Months Relative to the Administration of Treatment", side = 1, line = 1.385, at = 0.53,
      outer = TRUE, cex = 1.2)

mtext("Twenty Closest Control Units", side = 3, line = 1.385, at = 0.268,
      outer = TRUE, cex = 1.2)

mtext("Fifty Closest Control Units", side = 3, line = 1.385, at = 0.78,
      outer = TRUE, cex = 1.2)

mtext("Four Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.76,
      outer = TRUE, cex = 1.2)

mtext("Six Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.262,
      outer = TRUE, cex = 1.2)


dev.off()



###############################################################################################
## L = 4 ##
Maha_M10_ATE <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          log_transaction_area_l2 + 
                          log_transaction_area_l3 + 
                          log_transaction_area_l4 + 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 10,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M10_ATE, file = "output/Maha_M10_ATE.RData")

##
Maha_M5_ATE <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          log_transaction_area_l2 + 
                          log_transaction_area_l3 + 
                          log_transaction_area_l4 + 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 5,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M5_ATE, file = "output/Maha_M5_ATE.RData")

##
Pscore_Weighting_ATE <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 10,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_Weighting_ATE, file = "output/Pscore_Weighting_ATE.RData")


### Try L = 6 with Maha Matching
Maha_M10_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          log_transaction_area_l2 + 
                          log_transaction_area_l3 + 
                          log_transaction_area_l4 + 
                          log_transaction_area_l5 + 
                          log_transaction_area_l6 + 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 10,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M10_ATE_L6, file = "output/Maha_M10_ATE_L6.RData")

##
Maha_M5_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          log_transaction_area_l2 + 
                          log_transaction_area_l3 + 
                          log_transaction_area_l4 + 
                          log_transaction_area_l5 + 
                          log_transaction_area_l6 + 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Maha",
                        M = 5,
                        weighting = TRUE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Maha_M5_ATE_L6, file = "output/Maha_M5_ATE_L6.RData")

##
Pscore_M10_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 10,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M10_ATE_L6, file = "output/Pscore_M10_ATE_L6.RData")

Pscore_M5_ATE_L6 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 6, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 5,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -6:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M5_ATE_L6, file = "output/Pscore_M5_ATE_L6.RData")

load("output/Pscore_M5_ATE_L6.RData")

## 
Pscore_M10_ATE_L4 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 10,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M10_ATE_L4, file = "output/Pscore_M10_ATE_L4.RData")

Pscore_M5_ATE_L4 <- pforeach(i = 1:5, .cores = cores, .seed = 2018) ({
  matches <- PanelMatch(formula = log_transaction_area ~ 
                          tour + 
                          log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                          log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                          log_expd_l1 + log_middle_school_l1 + 
                          msec2currentsec_l1 + 
                          msec2currentgvn_l1 + 
                          mayor2currentsec_l1 + 
                          mayor2currentgvn_l1 + 
                          msec2currentsec2_l1_missing, 
                        lag = 4, max.lead = (i-1),
                        unit.id = "county_id",
                        time.id = "time",
                        treatment = "tour",
                        method = "Pscore",
                        M = 5,
                        weighting = FALSE,
                        qoi = "ate",
                        data = d_sub1)    
  results <- PanelEstimate(lead = -4:(i-1), matched_sets = matches,
                           qoi = "ate", inference = "bootstrap",
                           ITER = 500)
  out <- results
  list(out)
})

save(Pscore_M5_ATE_L4, file = "output/Pscore_M5_ATE_L4.RData")


########## Making Figures ##############
#### PS Figures ####

load("output/Pscore_M10_ATE_L4.RData")
load("output/Pscore_M5_ATE_L4.RData")
load("output/Pscore_M10_ATE_L6.RData")
load("output/Pscore_M5_ATE_L6.RData")

all_results <- list(Pscore_M10_ATE_L4, Pscore_M5_ATE_L4, Pscore_M10_ATE_L6, 
                    Pscore_M5_ATE_L6)


pdf(file = "output/Figure15.pdf", 
    width = 10, height = 8)
par(mfrow = c(2,2), 
    mar = c(2, 3, .5 , 1), oma = c(3, 3,4, 1))

for (i in 1:length(all_results)) {
  total_periods <- length(all_results[[i]][[length(all_results[[i]])]]$lead)
  P <- rep(NA, total_periods)
  L_90 <- rep(NA, total_periods)
  U_90 <- rep(NA, total_periods)
  L_95 <- rep(NA, total_periods)
  U_95 <- rep(NA, total_periods)
  CI_90 <- .9
  CI_95 <- .95
  if (i == 1|i == 3){
    ylab = "Estimated Effect of Inspection"
  } else {
    ylab = ""
  }
  for (j in 1:length(all_results[[i]])) {
    results <- all_results[[i]][[j]]
    if (j == 1){
      P[1:length(results$lead)] <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      L_90[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_90
      U_90[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      L_95[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_95
      U_95[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
    } else {
      tmp <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      P[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      U_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_95[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
      U_95[results$lag + j] <- tmp[length(tmp)]
    }   
  }
  grey_gov <- length(P) - 5 - .5       
  xlim_gov <- length(P) - .5
  P[length(P)-5] <- NA
  L_90[length(L_90)-5] <- NA
  U_90[length(U_90)-5] <- NA
  L_95[length(L_95)-5] <- NA
  U_95[length(U_95)-5] <- NA
  ylim <- c(min(L_95, na.rm = T) * 1.1, 
            max(U_95, na.rm = T) * 1.1) 
  plot(0:(length(as.numeric(P))-1), P, xlab = "", ylab = ylab, 
       xaxt = "n",
       # yaxt = "n",
       xlim = c(-.5,xlim_gov),
       
       panel.first = rect(-10, -10, grey_gov, 10, border = NA,
                          col = "grey88"), 
       ylim = ylim, 
       pch = 16, cex = 1.4) 
  #  }
  
  # #  if (k == 1) {
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_95[k+1], U_95[k+1]))
  }
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_90[k+1], U_90[k+1]), lwd = 5,
          col = "grey")
  }
  points(0:(length(as.numeric(P))-1), P, 
         pch = 16, cex = 1.4)
  abline(h = 0, lty = 3)
  axis(side = 1, at = 0:(length(P)-1), 
       labels = -(length(P)-5):4)
}

mtext("Months Relative to the Administration of Treatment", side = 1, line = 1.385, at = 0.53,
      outer = TRUE, cex = 1.2)

mtext("Five Closest Control Units", side = 3, line = 1.385, at = 0.268,
      outer = TRUE, cex = 1.2)

mtext("Ten Closest Control Units", side = 3, line = 1.385, at = 0.78,
      outer = TRUE, cex = 1.2)

mtext("Four Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.76,
      outer = TRUE, cex = 1.2)

mtext("Six Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.262,
      outer = TRUE, cex = 1.2)


dev.off()

##### Maha Figures #####
load("output/Maha_M10_ATE.RData")
load("output/Maha_M5_ATE.RData")
load("output/Maha_M10_ATE_L6.RData")
load("output/Maha_M5_ATE_L6.RData")

all_results <- list(Maha_M10_ATE, Maha_M5_ATE, Maha_M10_ATE_L6, 
                    Maha_M5_ATE_L6)


pdf(file = "output/Figure13.pdf", 
    width = 10, height = 8)
par(mfrow = c(2,2), 
    mar = c(2, 3, .5 , 1), oma = c(3, 3,4, 1))

for (i in 1:length(all_results)) {
  total_periods <- length(all_results[[i]][[length(all_results[[i]])]]$lead)
  P <- rep(NA, total_periods)
  L_90 <- rep(NA, total_periods)
  U_90 <- rep(NA, total_periods)
  L_95 <- rep(NA, total_periods)
  U_95 <- rep(NA, total_periods)
  CI_90 <- .9
  CI_95 <- .95
  if (i == 1|i == 3){
    ylab = "Estimated Effect of Inspection"
  } else {
    ylab = ""
  }
  for (j in 1:length(all_results[[i]])) {
    results <- all_results[[i]][[j]]
    if (j == 1){
      P[1:length(results$lead)] <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      L_90[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_90
      U_90[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      L_95[1:length(results$lead)] <-  t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                               results$o.coef, byrow = TRUE) - results$boots,
                                                      probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                      na.rm = T, drop = FALSE))[1,] # bc percentile CI_95
      U_95[1:length(results$lead)] <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                                              results$o.coef, byrow = TRUE) - results$boots,
                                                     probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                                                     na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
    } else {
      tmp <- 2*results$o.coef -
        colMeans(results$boots, na.rm = T)
      
      P[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_90)/2, CI_90+(1-CI_90)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_90
      
      U_90[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[1,] 
      
      L_95[results$lag + j] <- tmp[length(tmp)]
      
      tmp <- t(colQuantiles(2*matrix(nrow = results$ITER, ncol = length(results$o.coef), 
                                     results$o.coef, byrow = TRUE) - results$boots,
                            probs = c((1-CI_95)/2, CI_95+(1-CI_95)/2),
                            na.rm = T, drop = FALSE))[2,] # bc percentile CI_95
      
      U_95[results$lag + j] <- tmp[length(tmp)]
    }   
  }
  grey_gov <- length(P) - 5 - .5       
  xlim_gov <- length(P) - .5
  P[length(P)-5] <- NA
  L_90[length(L_90)-5] <- NA
  U_90[length(U_90)-5] <- NA
  L_95[length(L_95)-5] <- NA
  U_95[length(U_95)-5] <- NA
  ylim <- c(min(L_95, na.rm = T) * 1.1, 
            max(U_95, na.rm = T) * 1.1) 
  plot(0:(length(as.numeric(P))-1), P, xlab = "", ylab = ylab, 
       xaxt = "n",
       # yaxt = "n",
       xlim = c(-.5,xlim_gov),
       
       panel.first = rect(-10, -10, grey_gov, 10, border = NA,
                          col = "grey88"), 
       ylim = ylim, 
       pch = 16, cex = 1.4) 
  #  }
  
  # #  if (k == 1) {
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_95[k+1], U_95[k+1]))
  }
  for (k in 0:(length(as.numeric(P))-1)) {
    lines(c(k,k), c(L_90[k+1], U_90[k+1]), lwd = 5,
          col = "grey")
  }
  points(0:(length(as.numeric(P))-1), P, 
         pch = 16, cex = 1.4)
  abline(h = 0, lty = 3)
  axis(side = 1, at = 0:(length(P)-1), 
       labels = -(length(P)-5):4)
}

mtext("Months Relative to the Administration of Treatment", side = 1, line = 1.385, at = 0.53,
      outer = TRUE, cex = 1.2)

mtext("Five Closest Control Units", side = 3, line = 1.385, at = 0.268,
      outer = TRUE, cex = 1.2)

mtext("Ten Closest Control Units", side = 3, line = 1.385, at = 0.78,
      outer = TRUE, cex = 1.2)

mtext("Four Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.76,
      outer = TRUE, cex = 1.2)

mtext("Six Month Lags \n Estimated Effect of Inspection", side = 2, line = -.5, at = 0.262,
      outer = TRUE, cex = 1.2)

dev.off()


